Prevalence and distribution of Varroa destructor and Nosema spp. in symptomatic honey bee colonies across the USA from 2015 to 2022

USDA-ARS Bee Research Laboratory received symptomatic honey bee (Apis mellifera L.) samples across the United States for disease diagnosis. Here, we present a retrospective study and cartography of ectoparasite Varroa destructor and intracellular microsporidia parasite Nosema spp. These two major parasites were identified in the diseased honey bee samples between 2015 and 2022. Varroa infestation level (VIL) was examined by a wash technique (Mites/100 bees) and calculated as a percentage, while Nosema infection was quantified by microscopical spore count (Million Spores/Bee). Data were analyzed by month, year, state, and by nine geographical climate regions described in the U.S. Of adult bee samples (n = 4039) that were analyzed for Varroa mite infestation, the overall VIL in the U.S. ranged between 0.4 and 30.85%, with an overall national VIL and Varroa prevalence of 8.21% and 85.14%, respectively. Overall monthly data showed VIL constantly exceeded the critical level of 4% except from June to September and reached a maximum of 15% in January and December. Nationwide, VIL significantly (p < 0.001) increased from 2015 to 2018 (1.1–4.7%), plateaued from 2018 to 2021 (4.7–4.5%), followed by a significant decrease in 2022 (3.6%). Significant VIL differences (p < 0.001) were recorded among climate regions, with the highest mite infestation levels in the Upper Midwest region (13.9%) and the lowest in the West region (5.1%). Of adult bee samples (n = 2,994) that were analyzed for Nosema infection, Nosema spore count ranged between (1–16.8) million spores per bee among states, with a national average of 6.8 and a prevalence of 99.7%. The lowest and highest Nosema loads were respectively recorded in the South region (3.1) and Upper Midwest (10.5), a significant difference (p < 0.001). No statistical differences were recorded among the six other climate regions. Overall, VIL and Nosema infection correlated significantly (p < 0.001) with a regression coefficient of (R2 = 0.6). Our data, which originated from ailing bee colonies, showed significantly higher rates of maladies compared to data from healthy colonies obtained by the USDA-APHIS National Honey Bee Survey, demonstrating the role of bee diseases caused by Varroa mite and Nosema in honey bee population declines.


VIL and climate regions
Our results evidenced statistically significant (p < 0.001) variations in Varroa infestations among climate regions, whether per month or as overall averages (Fig. 5).In the monthly load analysis, the lowest elbow of the VIL curves matches the period between June to August.The highest overall VIL was recorded in the Upper Midwest region (IA, MI, MN, WI), while the lowest VIL was documented in the West region (CA and NV) (Fig. 5a and Table 1).VIL of other climate regions ranged between these two ends, with significant inter-region differences detailed in Fig. 5a.Concerning the NHBS data, monthly lower VILs per region are recorded between March to July, with significant increases in Autumn (Sep-Nov) (Fig. 5b).The highest overall VIL in the NHBS data was recorded in the Ohio Valley region, while the lowest VIL was in the West region, similar to what was obtained in our data (Fig. 5b).
Vol:.( 1234567890) Analysis of Nosema spore load by month from different climate regions showed significant differences (p < 0.001) among the latter (Fig. 10a).Highest levels of infection were recorded from January to April, followed by declines from May to November and increased again in December (Fig. 10a).The highest overall Nosema spore count was recorded in Upper Midwest region and the lowest in the South region.No significant differences in Nosema load were found among the other six climate regions (Fig. 10b).

VIL and Nosema relationship
A regression analysis was conducted to elaborate on the relationship between mite load and Nosema spore count.This analysis was carried out on national monthly averages for both variables.It showed a significant (p < 0.001) positive linear trend (R 2 = 0.6) between VIL and Nosema infection in our analyzed samples (Fig. 10c).Similarly, multiple positive yearly correlations between VIL and Nosema infection were identified, as well as two negative correlations between Nosema loads in 2021 and both 2016 and 2019 with R = − 0.7 and − 0.6 respectively (Fig. 10d).

Discussion
The Varroa destructor, considered one of the most devastating parasites and the leading cause of colony mortality, was identified in all state samples analyzed over eight years.Few previous studies from Latin American countries reported VIL-related data.In Colombia, mainly among Africanized honey bee (AHB) populations, a VIL of 4.5% was reported 34 (Table 2).Mexico, which comprises a mix of Africanized honey bees (AHB) and European honey bees (EHB) populations 35,36 , reported VIL ranging between 3.5 and 7.4% 37,38 , Table 2.A recent longitudinal study from Brazil reviewing VIL from 1979 to 2020 within AHB populations reported a VIL of 4.5% 39 .Similarly close VIL of 5.4% was reported in the same country in another study 40 .The U.S. national VIL found within clinically diseased colonies (8.21%) fits closer to the end range found in Mexico, which relies on both EHB and AHB populations (Table 2).Understanding the genetic background of the honey bee populations is essential and critical in this context, as AHB populations were described to be more resistant to Varroa mite infestation compared to EHB populations [41][42][43] .Arguably one of the most comprehensive U.S. data related to mite infestation and Nosema spore load was generated by the USDA-APHIS National Honey Bee Survey (NHBS).This survey has been conducted on healthy honey bee colonies across the U.S.A. since 2009.In order to explore further our findings, we compared our data with those of the NHBS within the same year range (2015-2020) of our current study generated at the Bee Research Laboratory (BRL).Clearly, the samples we received and analyzed reflect clinically diseased populations or colonies as both overall VIL and Nosema loads are significantly higher than what was found in the data from the NHBS, which was randomly conducted nationwide on healthy colonies (Figs.S1, S2).
Interestingly, Varroa prevalence was higher (91.52%) in the NHBS data than our BRL data (85.14%),indicating significant proliferation of this parasite regardless of the colony's health status (Fig. S1b).This was not the case for Nosema prevalence, which was approximately two-fold higher in our samples (99.73%) than in healthy colonies (46.54%) sampled by the NHBS (Fig. S2b).The overall VIL recorded in our diagnosed samples (8.21%) was significantly higher than the nationwide average VIL of healthy colonies reported by the NHBS (3.09%) (Fig. S1).Therefore, it is conceivable to hypothesize that elevated levels of Varroa mites may have triggered the decline of these colonies.Moreover, monthly Varroa load in diseased, weak, or struggling colonies did not follow the same patterns identified in healthy colonies (Fig. 4).For instance, samples from weak and/or dying colonies shipped to us in winter time had the highest VIL.
In contrast, the NHBS data, which displays the classic bee-Varroa parasitism dynamic, showed a single high peak in October to drop again during winter.Hence, the monthly VIL pattern of the diseased samples (BRL data) could reflect the intensification of Varroa treatments in unwell colonies during the beekeeping season (Apr-Oct), which had pushed the parasite development and survival to Winter, when beekeepers usually apply no treatment.In such struggling colonies, a take-over of the parasite has likely been established, leading to bees' failure to suppress Varroa mites.Such a phenomenon is supported by a constantly elevated VIL exceeding the 4% critical threshold 44 in our samples.On the other hand, the host-parasite dynamic in healthy colonies differed from weak declining colonies as the former's social resilience and ability to suppress diseases remained intact.This is supported by the fact that colonies of the NHBS maintained a relatively low range of VIL (1-6%) across the year (Fig. 4).Eventually, such resilience towards Varroa mites was not due to natural resistance per se, but rather to intensive pest management control carried out by beekeepers nationwide.BRL and NHBS data firmly agreed that colonies in the West region (CA and NV) had the lowest VIL.Despite slight differences between both datasets vis-à-vis regions with the highest VIL, the consensus indicates that northern zones with colder climates favored higher VIL.Besides significant variations in VIL among climate regions, the absence of meaningful correlations of VIL among regions is strong evidence of the implication of climate region on Varroa mite infestation.The highest Nosema infestation identified in the Upper Midwest region reveals a more effective proliferation and intensification of this pathogen in cold, moist environments, as pointed out in a previous study 45 .1.
In conclusion, this study provides valuable information and significant insights into the seasonal variation and distribution of Varroa mite infestation and Nosema infection among diseased honey bee colonies in the U.S.A.This underscores the necessity for ongoing commitments from beekeepers and researchers to develop and implement strategies to prevent and mitigate honey bee diseases and stressors, including integrated pest management practices, breeding for disease-resistant bees, and promoting habitat conservation.Addressing honey bee disease challenges is of utmost significance to safeguard the well-being of honey bees and to uphold their pivotal role in providing essential pollination services.

Honey bee samples
The Bee Research Laboratory (BRL) received a total of 7,033 symptomatic honey bee samples between 2015 and 2022.Honey bee samples are usually soaked in 70% ethanol and shipped via post mail.Based on our guidelines, State Apiary Inspectors and beekeepers typically collected approximately 300 adult bees from each colony with signs of diseases, weakening or declining colonies, and in some cases, dead bees from collapsed colonies.The identification information, including the sampling date, locations (state, zip code), and names of beekeepers or apiary inspectors, were all marked on the plastic zipper bags used for holding bee samples.

Varroa infestation level
The Varroa Infestation Level (VIL) was described for each sample based on the number of mites identified.VIL was calculated as a percentage by dividing the number of mites by the number of bees in a sample and multiplying by 100 as described in previous studies 34,46 .Varroa mite load per bee sample was obtained using a multistep process.Firstly, the bees in each sample were agitated and transferred into a screening container fitted with a #8 mesh insert.This mesh effectively retained the bees while allowing debris and varroa mites to pass through.The screening container was placed within a fine mesh sieve and placed underneath a tap.The tap was opened fully to pressure wash the bees for 40 s while continuously swirling the screening container.Once the washing process was completed, the screening container was set aside for water to drain off, and any remaining water  www.nature.com/scientificreports/ was absorbed using tissue paper.The total number of bees per sample was determined based on the weight of a single wet bee (0.16129 g), which was obtained from previous trials and regression correlations conducted in our lab.Subsequently, each sample's total number of bees was calculated based on their collective weight.The mites collected on the fine mesh during the washing process were carefully counted, and VIL was calculated as described above.

Nosema infection
For each colony, 30 honey bee workers were randomly subsampled and examined.Honey bee abdomens were removed and placed in a zip-lock bag.The abdomens were crushed within the bag using a rolling pin on a flat surface.Subsequently, 30 mL of distilled water was added to the bag to create a homogenous suspension.Microscopic examination was conducted on 10 µL of the suspension, which was transferred to a hemocytometer, examined under a light microscope at 400× magnification 47 , and screened for the presence or absence of Nosema spp.If Nosema spores were observed, they were quantified, and the infestation rates of Nosema spp.were calculated based on the number of spores counted per examined field 48 , with the results expressed as millions of spores per bee.

VIL and Nosema spore loads across the U.S. climate regions
This prevalence of Varroa and Nosema is expressed as the percentage of samples that tested positive for either of these parasites within a specific state or geographic region.Varroa and Nosema prevalence were analyzed and studied from a geographical perspective based on climate regions.The classification of climate regions used in this study was based on the official U.S. government classification, as outlined by the National Centers for Environmental Information.According to this classification, the U.S. states were grouped into nine major climate regions 49 , namely: (1) Northeast (NE), (2) Northern Rockies and Plains (NRP), (3) Northwest (NW), (4) Ohio Valley (OV), ( 5) Southeast (SE), ( 6) Southwest (SW), (7) Upper Midwest (UMW) and ( 8) West (W).Further information regarding the specific states within each climate region is provided in Table 1, along with the corresponding abbreviations used throughout the text and figures.

National honey bee survey NHBS data
Since 2009, the USDA-APHIS, through the National Honey Bee Survey program (NHBS), has conducted an ongoing random and regular honey bee disease screening on healthy colonies across the U.S. As a means of

Data analysis
Statistical analyses and statistically-related figures for each dataset were carried out in the R environment 50 via RStudio 51 Version (2022.07.0).Geographical mapping of honey bee pests was conducted using the Tableau Public platform (https:// www.table au.com).All datasets were tested for normality using the Shapiro test before conducting statistical analyses.Kruskal-Wallis rank test, a non-parametric test, was conducted at a 95% confidential interval with three levels of significance (p < 0.05*; < 0.001**; < 0.001***) on data that failed the normality test.Adjustments to p-values were made where appropriate using the Benjamini-Hochberg method.Figures were generated in the R environment utilizing three main Libraries: "ggplot2", "doby", and "plyr".Mite and Nosema loads were displayed longitudinally per month and as overall averages per year or climate region.Error bars of all line graphs represent the Standard Error (SE) except for the boxplots (box and whisker plots), which display the median, first and third quartiles, and both maximum and minimum values of variables.Outliers from Box plot graphs were omitted in some instances with log transformation for better data visualization.Average VIL and Nosema spore count per state and nationwide, summarized mainly in Tables, represent the mean values and not the medians unless stated otherwise.Heatmaps were generated for VIL data using the "pheatmap" library to visualize the relationships between mite load vis-à-vis month and climate region.Correlation analyses between variables were all conducted at a cutoff of p < 0.05 with the "spearman" method as the data failed normality utilizing both the "corrplot" and "Hmisc" Libraries.Linear regression analysis was conducted to study the relationships between both Varroa and Nosema loads as variables.

Figure 1 .
Figure 1.Overall number and distribution of analyzed samples from 2015 to 2022 across the U.S., displayed per state and county Zip code.(a) Samples analyzed for Varroa mite per 100 bees and (b) Nosema spore count per bee.Overlapping samples in the Zip code maps were omitted for better visualization of the sample distribution.

Figure 3 .
Figure 3. Percentage of average Varroa Infestation Level (VIL), displayed per state and year (2019-2022).States with gray font did not send samples for analysis.Mites were counted from adult honey bees using a wash technique as described in the materials and methods.VIL range and Varroa prevalence (VP) per year are summarized on top of each map.

Figure 4 .
Figure 4. Nationwide averages of Varroa Infestation Level (VIL) displayed by month and year for data of the current study (a) and the USDA-APHIS National Honey Bee Survey NHBS (b).(n) is the number of samples in each graph and statistical analysis.Error bars in the line graphs are the Standard Error (SE).The Kruskal-Wallis test was conducted at three levels of significance: p < 0.05*, p < 0.01**, p < 0.001***.The box plot's median is displayed in white font, and box plots with different alphabetic letters are statistically significant.

Figure 5 .
Figure 5. Nationwide Varroa Infestation Level (VIL) displayed by climate region per month and overall average on data of the current study (a) and the USDA-APHIS National Honey Bee Survey NHBS (b).(n) is the number of samples in each graph and statistical analysis.Error bars in the line graphs are the Standard Error (SE).The Kruskal-Wallis test was conducted at three levels of significance: p < 0.05*, p < 0.01**, p < 0.001***.Box plots with different alphabetic letters are statistically significant.Climate region abbreviations are given in Table1.

Figure 6 .Figure 7 .
Figure 6.Heatmap and correlation Varroa Infestation level (VIL) per month and climate region.VIL heatmap per month and climate region and correlation of VIL among climate regions on current study data (a) and the NHBS data (b).White crosses are unavailable values from NHBS data.Correlation analyses were conducted at a cutoff of p < 0.05, R-values are given within pairwise correlation circles, only significant correlations are displayed, and blank squares were not statistically significant.The assignment of states to climate regions is detailed in Table1.

Figure 8 .
Figure 8. Overall averages of Nosema spores identified per bee (Million Spores/Bee), displayed per state and year (2019-2022).States with gray font did not send samples for analysis.Nosema spore counts were conducted on adult honey bees as described in the materials and methods.National Nosema range and prevalence (NP) per year are summarized on top of each map.

Figure 9 .
Figure 9. Nationwide Nosema spore count displayed by month and year for both datasets: Current study (a) and the National Honey Bee Survey NHBS (b).(n) is the number of samples in each graph.Error bars in the line graphs are the Standard Error (SE).The Kruskal-Wallis test was conducted at three levels of significance: p < 0.05*, p < 0.01**, p < 0.001***.The box plot's median is displayed in white font, and box plots with different alphabetic letters are statistically significant.

Figure 10 .
Figure 10.Climate region implication on Nosema load and correlations between VIL and Nosema infection.(a) Monthly longitudinal display of average spores per climate region.(b) Overall spore average per climate region.(c) Regression line and coefficient of Varroa mite load and Nosema spore count.(d) Correlation of Varroa mite infestation and Nosema infection conducted by average year load.The correlation was conducted at p < 0.05 cutoff, R-values are given within pairwise correlation circles, only significant correlations are displayed, and blank squares were not significant.Box plots with different alphabetic letters are statistically significant.

Table 1 .
Overall nationwide and state summary of Varroa mite infestation and Nosema spore count between 2015 and 2022.N1: Number of analyzed samples per state for Varroa mite, N2: Samples analyzed for Nosema.T1: Total number of samples analyzed per state.T2: National sums and averages, VP: Varroa prevalence, NS: Nosema spore count (Million Spores/Bee), and NP: Nosema prevalence.States were assigned to climate regions according to the U.S. National Centers for Environmental Information (Karl and Koss 1984).